// ///////////////////////////// MIT License //////////////////////////////////// //
//                                                                                //
// Copyright (c) 2010 David Zsolt Manrique                                        //
//                    david.zsolt.manrique@gmail.com                              //
//                                                                                //
// Permission is hereby granted, free of charge, to any person obtaining a copy   //
// of this software and associated documentation files (the "Software"), to deal  //
// in the Software without restriction, including without limitation the rights   //
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      //
// copies of the Software, and to permit persons to whom the Software is          //
// furnished to do so, subject to the following conditions:                       //
//                                                                                //
// The above copyright notice and this permission notice shall be included in     //
// all copies or substantial portions of the Software.                            //
//                                                                                //
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     //
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         //
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  //
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN      //
// THE SOFTWARE.                                                                  //
//                                                                                //
// ////////////////////////////////////////////////////////////////////////////// //

#include <iostream>
#include <ostream>
#include <sstream>
#include <vector>
#include <algorithm>
#include <complex>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <fstream>
#include <string>
#include <algorithm>
#include <iomanip>
#include <vector>
#include <list>
#include <set>
#include <cstdarg>
#include <stdexcept>

#include <Eigen/Eigen>

#include <cmdioutil.h>


void get_atomic_block_IJ(const Eigen::MatrixXcd & full,const std::vector<std::pair<int,int> > & orbital_range,int I,int J,Eigen::MatrixXcd & block)
{
    int i,j,ni,nj;
    i = orbital_range[I].first;
    ni = orbital_range[I].second;
    j = orbital_range[J].first;
    nj = orbital_range[J].second;
    block = full.block(i,j,ni,nj);
}

void set_atomic_block_IJ(Eigen::MatrixXcd full,const std::vector<std::pair<int,int> > & orbital_range,int I,int J,const Eigen::MatrixXcd & block)
{
    int i,j,ni,nj;
    i = orbital_range[I].first;
    ni = orbital_range[I].second;
    j = orbital_range[J].first;
    nj = orbital_range[J].second;
    full.block(i,j,ni,nj) = block;
}

void get_block(const Eigen::MatrixXcd & full,const std::vector<std::pair<int,int> > & orbital_range,const std::list<int> & index_list_I,const std::list<int> & index_list_J,Eigen::MatrixXcd & block)
{
    int nI,nJ;
    nI = 0;
    for(std::list<int>::const_iterator it = index_list_I.begin(); it != index_list_I.end(); it++) nI += it->second;
    nJ = 0;
    for(std::list<int>::const_iterator it = index_list_J.begin(); it != index_list_J.end(); it++) nJ += it->second;
    
    block.resize(nI,nJ);
    int li = 0;
    for(std::list<int>::const_iterator it = index_list_I.begin(); it != index_list_I.end(); it++) 
    {
        int I = *it;
        int i,ni;
        i = orbital_range[I].first;
        ni = orbital_range[I].second;
        int lj = 0;
        for(std::list<int>::const_iterator jt = index_list_J.begin(); jt != index_list_J.end(); jt++) 
        {
            int J = *jt;
            int j,nj;
            j = orbital_range[J].first;
            nj = orbital_range[J].second;
            block.block(li,lj,ni,nj) = full.block(i,j,ni,nj);
            lj += nj;
        }
        li += ni;
    }
}


void set_block(const Eigen::MatrixXcd & full,const std::vector<std::pair<int,int> > & orbital_range,const std::list<int> & index_list_I,const std::list<int> & index_list_J,Eigen::MatrixXcd & block)
{
    int li = 0;
    for(std::list<int>::const_iterator it = index_list_I.begin(); it != index_list_I.end(); it++) 
    {
        int I = *it;
        int i,ni;
        i = orbital_range[I].first;
        ni = orbital_range[I].second;
        int lj = 0;
        for(std::list<int>::const_iterator jt = index_list_J.begin(); jt != index_list_J.end(); jt++) 
        {
            int J = *jt;
            int j,nj;
            j = orbital_range[J].first;
            nj = orbital_range[J].second;
            full.block(i,j,ni,nj) = block.block(li,lj,ni,nj)
            lj += nj;
        }
        li += ni;
    }
}


void load_hamiltonian(std::string fn,std::vector<std::pair<int,int> > & range,
                                    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> & h,
                                    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> & s,double &ef_m)
{
    std::complex<double> ii(0,1.0);
    
    std::vector<std::string> lines;
    io::load(fn,"#",lines);

    ef_m = 0.0;
    std::vector<std::string> efermi_block;
    io::parse_block("<fermi-energy>","</fermi-energy>",lines,efermi_block);
    for(std::vector<std::string>::iterator it = efermi_block.begin(); it != efermi_block.end(); it++)
    {
        std::istringstream iss(*it);
        if (!(iss >> ef_m)) throw std::runtime_error("Not able to parse fermi energy block in the input file!");
        break;
    }

    std::vector<std::string> orbital_block;
    if(!io::parse_block("<index-orbital>","</index-orbital>",lines,orbital_block)) throw std::runtime_error("No or inproper index-orbital block in the input file!");
    if(!orbital_block.size()) throw std::runtime_error("No or inproper index-orbital block in the input file!");
    
    range.clear();
    int index,orbital;
    int inc = 0;
    int orbital_index = 0;
    for(std::vector<std::string>::iterator it = orbital_block.begin(); it != orbital_block.end(); it++)
    if(util::match("% %",*it,index,orbital))
    {
            std::pair<int,int> p(orbital_index,orbital);
            orbital_index += orbital;
            range.push_back(p);
            if(++inc != index) throw std::runtime_error("Index has to be in order from 1 in the index-orbital block!");
    }
    int n = orbital_index;
    h.resize(n,n); h.setZero(); s = h; s.setIdentity();
    double re,im;
    int i,j;

    std::vector<std::string> overlap_block;
    if(!io::parse_block("<overlap>","</overlap>",lines,overlap_block)) throw std::runtime_error("There is no overlap matrix block or it is inproper in the input file!");
    for(std::vector<std::string>::iterator it = overlap_block.begin(); it != overlap_block.end(); it++)
    if(util::match("% % % % ",*it,i,j,re,im)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) s(i-1,j-1)=re+ii*im; else throw std::runtime_error("Inconsistent overlap matrix index with index orbital block!");  }
    else if(util::match("% % % ",*it,i,j,re)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) s(i-1,j-1)=re+ii*0.0; else throw std::runtime_error("Inconsistent overlap matrix index with index orbital block!"); }
    
    std::vector<std::string> hamiltonian_block;
    if(!io::parse_block("<hamiltonian>","</hamiltonian>",lines,hamiltonian_block)) throw std::runtime_error("There is no hamiltonian block or it is inproper in the input file!");
    for(std::vector<std::string>::iterator it = hamiltonian_block.begin(); it != hamiltonian_block.end(); it++)
    if(util::match("% % % % ",*it,i,j,re,im)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) h(i-1,j-1)=re+ii*im; else throw std::runtime_error("Inconsistent hamiltonian index with index orbital block!"); }
    else if(util::match("% % % ",*it,i,j,re)) { if((i<n+1)&&(j<n+1)&&(i>0)&&(j>0)) h(i-1,j-1)=re+ii*0.0; else throw std::runtime_error("Inconsistent hamiltonian index with index orbital block!"); }
    if(!hamiltonian_block.size()) throw std::runtime_error("There is no hamiltonian block or it is inproper in the input file!");

}




void load_hamiltonian_gollum(std::string fn,std::vector<std::pair<int,int> > & range,
                                    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> & h,
                                    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> & s,double &ef_m)
{
    std::ifstream is(fn.c_str());
    std::complex<double> ii(0,1.0);
    std::string dump;
    is >> dump >> dump >> dump >> dump >> dump >> dump >> dump;
    is >> dump >> dump >> dump >> dump >> dump >> dump >> ef_m;
    is >> dump >> dump >> dump >> dump >> dump >> dump;
    int rows,cols; 
    is >> dump >> dump >> rows;
    is >> dump >> dump >> cols;
    int index;
    int pindex = 1;
    int pi = 0;
    int n = rows;    
    for(int j = 0; j < n; j++)
    {
        is >> index >> dump >> dump >> dump >> dump;
        if(index != pindex)
        {
            std::pair<int,int> p(pi,j-pi);
            range.push_back(p);
            pindex = index;
            pi = j;
        }
    }
    std::pair<int,int> p(pi,n-pi);
    range.push_back(p);

    is >> dump >> dump >> dump >> dump >> dump >> dump;
    //int rows,cols;
    is >> dump >> dump >> rows;
    is >> dump >> dump >> cols;
    for(int j = 0; j < rows; j++)
        is >> dump >> dump >> dump;

    is >> dump >> dump >> dump >> dump >> dump >> dump;
    //int rows,cols;
    is >> dump >> dump >> rows;
    is >> dump >> dump >> cols;
    h.resize(n,n); h.setZero(); s = h;
    for(int k = 0; k < rows; k++)
    {
        int i,j;
        double res,ims,reh,imh;
        is >> dump >> i >> j >> res >> ims >> reh >> imh;
        h(i-1,j-1)=reh+ii*imh;
        s(i-1,j-1)=res+ii*ims;

    }
    is.close();
    
}


#ifdef DEBUG
void out_complex_matrix(Eigen::Matrix< std::complex<double> , Eigen::Dynamic, Eigen::Dynamic >  m,std::ostream & os,bool comment = true,int w = 18,int p = 10)
{
        //os.flags(std::ios::scientific);
        os.flags(std::ios::fixed);
        if(comment)
        {
                os << '#';
                os << std::left << std::setw(w-1) << std::setprecision(p) << "Re Im ... Re Im" << ' '; 
                os << std::endl;
        }
        for(int i = 0; i < m.rows(); i++)
        {
	        for(int j = 0; j < m.cols(); j++)
                {
                        os << std::right << std::setw(w+4) << std::setprecision(p) << real(m(i,j)) << ' '; 
                        os << std::right << std::setw(w) << std::setprecision(p) << imag(m(i,j)) << ' '; 
                }
                os << std::endl;                        
        }                
        os << std::endl;                        
}
#endif

#ifdef DEBUG
void out_complex_vector(Eigen::Matrix< std::complex<double> , Eigen::Dynamic, 1 >  m,std::ostream & os,bool comment = true,int w = 18,int p = 10)
{
        os.flags(std::ios::scientific);
        if(comment)
        {
                os << '#';
                os << std::left << std::setw(w-1) << std::setprecision(p) << "Re Im" << ' '; 
                os << std::endl;
        }
        for(int i = 0; i < m.rows(); i++)
        {
                        os << std::right << std::setw(w+4) << std::setprecision(p) << real(m(i)) << ' '; 
                        os << std::right << std::setw(w) << std::setprecision(p) << imag(m(i)) << ' '; 
                os << std::endl;                        
        }                
        os << std::endl;                        
}
#endif

#ifdef DEBUG
void out_real_matrix(Eigen::Matrix< double , Eigen::Dynamic, Eigen::Dynamic >  m,std::ostream & os,bool comment = true,int w = 18,int p = 10)
{
        os.flags(std::ios::scientific);
        for(int i = 0; i < m.rows(); i++)
        {
	        for(int j = 0; j < m.cols(); j++)
                {
                        os << std::right << std::setw(w+4) << std::setprecision(p) << m(i,j) << ' '; 
                }
                os << std::endl;                        
        }                
        os << std::endl;                        
}
#endif



void get_orbital_range(const std::vector<std::pair<int,int> > & range,const util::pair_expression<int,int,':'> &expr, int &io,int &no)
{
        int index1,index2,n;
        n = range.size();
        index1 = std::get<0>(expr);
        index2 = std::get<1>(expr);
        if(index1*index2 == 0) throw std::runtime_error("Index cannot be zero.");
        int i1,i2;
        i1 = index1 < 0 ? index1 + n : index1 - 1;
        i2 = index2 < 0 ? index2 + n : index2 - 1;
        if(i1 > i2) std::swap(i1,i2);                    
    
        io = range[i1].first;
        no = range[i2].second-io;
}


void structure_hkl(const std::vector<Eigen::Vector3d> & coord,const Eigen::Vector3d & a1,const Eigen::Vector3d & a2,const Eigen::Vector3d & a3,Eigen::MatrixXi & h,Eigen::MatrixXi & k,Eigen::MatrixXi & l)
{
    Eigen::Matrix3d ic,c;
    c.col(0) = a1;
    c.col(1) = a2;
    c.col(2) = a3;
    
    ic = c.inverse();
    
    int n = coord.size();
    h.resize(n,n);
    k.resize(n,n);
    l.resize(n,n);
    
    Eigen::Vector3d hkl;
    
    for(int i=0;i<n;i++)
    for(int j=0;j<n;j++)
    {
       hkl=ic*(coord[i]-coord[j]);
       h(i,j)=std::round(hkl.x());    
       k(i,j)=std::round(hkl.y());    
       l(i,j)=std::round(hkl.z());    
    }
}

void get_cmd(int argc, char* argv[],std::vector<Eigen::Vector3d> & coord,std::vector<std::pair<int,int> > & range,Eigen::VectorXcd & hfull,Eigen::VectorXcd & sfull,double & ef,std::vector<Eigen::Vector3d> & coord_model)
{
        cmd::arguments arg(argc,argv);
        cmd::switcharg<false> help(arg,"-h","--help", "print help message");
        cmd::switcharg<false> info(arg,"-i","--info", "print info/example message");

        cmd::switcharg<false> gollum(arg,"-g","--gollum", "hamiltonian file in gollum format");

        cmd::posarg<std::string> xyz_fn(arg,"reference structure input file name", "");
        cmd::posarg<std::string> ham_fn(arg,"reference hamiltonain input file name", "");
        cmd::posarg<std::string> model_xyz_fn(arg,"model structure input file name", "");

        if(*help)
        {
                arg.print_help(std::cout,"hamanal");
                std::cout << "examples:" << std::endl;
                std::cout << "version 0.9b" << std::endl;
                exit(EXIT_SUCCESS);
        }
        if(*info)
        {
                std::cout << "version 0.9b" << std::endl;
                std::cout << "Copyright (c) 2010 David Zsolt Manrique (david.zsolt.manrique@gmail.com)"  << std::endl;
                //return EXIT_SUCCESS;
                exit(EXIT_SUCCESS);
        }                                                                                                                
        
        if(!~xyz_fn) throw std::runtime_error("Give reference structure file!");
        if(!~ham_fn) throw std::runtime_error("Give reference hamiltonian file!");
        if(!~model_xyz_fn) throw std::runtime_error("Give model structure file!");
        
        xyz c1,c2;
        load(*xyz_fn,c1);
        load(*model_xyz_fn,c2);
        
        ef=0.0;
        if(!*gollum) load_hamiltonian(*ham_fn,range,hfull,sfull,ef);
        else load_hamiltonian_gollum(*ham_fn,range,hfull,sfull,ef);

}


